Sharp transition for single polarons in the one-dimensional Su-Schrieffer-Heeger model 
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We study a single polaron in the Su-Schrieffer-Heeger (SSH) model using four different techniques 
(three numerical and one analytical). Polarons show a smooth crossover from weak to strong cou- 
pling, as a function of the electron-phonon coupling strength A, in all models where this coupling 
depends only on phonon momentum q. In the SSH model the coupling also depends on the electron 
momentum k; we find it has a sharp transition, at a critical coupling strength A c , between states 
with zero and nonzero momentum of the ground state. All other properties of the polaron are also 
singular at A = A c . This result is representative of all polarons with coupling depending on k and q, 
and will have important experimental consequences (eg., in ARPES and conductivity experiments). 
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PACS numbers: 72.10.-d, 71.10.Fd, 71.38.-k 

Polarons have been of broad interest in physics ever 
since they were introduced in 1933, to describe dielec- 
tric charge carriers [l[. Apart from their central role in 
solid-state physics, with many models now in use 0d3|; 
they exemplify in quantum field theory the passage from 
weak to strong coupling in a non-trivial model of a single 
particle coupled to a bosonic field 0. The first serious 
non-perturbative studies by Feynman 0] of the Frohlich 
polaron, are now a classic, but only recently were accu- 
rate results established across the whole range of cou- 
pling strengths [7j. Since then, exact numerical studies 
have been made of, eg. , D-dimensional Holstein polarons 
in various lattice geometries, with D — 1,2,3 [8[; of 3D 
Rashba-Pekar polarons with short-range interactions [t| ; 
of pseudo Jahn- Teller polarons [l(| ; and so on. 

A central question in this field has been whether a 
sharp transition can exist in the polaronic ground state 
as a function of the dimensionless effective particle-boson 
coupling A. In all the above-cited work there is simply a 
smooth crossover, expected when the coupling depends 
only on the bosonic momentum q; then there must always 
be non-zero matrix elements between the ground state 
and excited polaron states ll|. However, quite generally, 
one expects the coupling to depend on both q and the 
particle momentum k; and then much less is known. 

In this paper we study a specific example of this general 
case. The particle-boson coupling is taken from the well- 
known "SSH model", introduced to describe electrons 
in 1-d polyacetylene 12j. Here we focus on the single 



polaron limit, not the more common case of half- filling, 
and the bosons are chosen to describe optical phonons. 
While this ignores the acoustic phonons which exist in 
real materials, it allows a direct comparison with the 
large number of results known for models which have 
a purely g-dependent coupling. The Hamiltonian thus 



takes the simple form % = H Q + V + H p h , where 



H = -t Q y^(c|c i+ x + h.c.) = 



phi 

^e k c\c k , (1) 
k 



describes the hopping of electrons between sites, with 
band dispersion e k = — 2i D cos(fc) (cj creates an electron 
on site i; cl creates a momentum state k). The term 
H p h — uiph b\bi describes dispersionless phonons {b\ 
creates a phonon on site i). The interaction is 



V 



x i+1 )(c\a + i 



h.c.) 



= N- i ^j2 M (k,M + M b -< ! + b <i) 



(2) 



k,q 



with site displacements Xi = 2 Mui b ij > anc ^ an 

interaction vertex 

M(k,q) = 2ia[sin(fc + q) - sin(fc)] 

= i(2\Lu ph t ) 1/2 [shi(k + q) - sin(fc)] (3) 

This interaction, with associated energy a = 
cW°y5ig^rT' describes the modulation of the hop- 
ping amplitude by phonons. We henceforth set 
t = 1, and define two dimensionless param- 
eters: the electron-phonon coupling parameter 
A = 2a 2 /{t o L0p h ) = (\M(k,q)\ 2 )/(2t Lu ph ), where (•) 
averages over the Brillouin zone; and the 'adiabaticity' 
ratio ujph/to (= u p h when t = 1). 

(i) Results: We treat this non-perturbative problem 
with the Momentum Average (MA) analytical approx- 
imation [(J [TH [l4| and three different numerical tech- 
niques: the Diagrammatic Monte Carlo (DMC) pj, the 
Limited Phonon Basis Exact Diagonalization (LPBED) 
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FIG. 1: (Color online) The polaron dispersion relation E(k) — 
E(k — 0) is shown in (a) and (c), and the GS Z-factors Z(k) 
at momentum k are shown in (b) and (d). Red (blue) tri- 
angles (circles) correspond to LPBED (BDMC) methods. In 
(a), (b), where uj ph = 0.5, A = 0.25,0.5,1.0,1.094,1.21,1.96 
(from top to bottom. In (c), (d), where ui p h = 3, A = 
0.25,0.5,1.0,2.0,4.0 (from top to bottom). MA results are 
shown as green solid curves. In (e) the GS energy for 
ijJph = 0.5 (upper line) and uj p h = 3 (lower curve) is shown; 
triangles, rhombi, squares and circles correspond to LPBED, 
MA, DMC, and BDMC methods, respectively. 

161 , and the Bold Diagrammatic Monte Carlo (BDMC) 
17] methods. Applications of the first three methods 
to polaron problems are well documented. However, our 
implementation of the BDMC method for the SSH model 
contains several new elements, reviewed in the support- 
ing material. 

In the following we display results as functions of A in 
both the adiabatic regime (choosing ui P h = 0.5), and the 
non-adiabatic regime (choosing u: p h — 3.0). We begin 
with the quasiparticle dispersion E(k) and renormaliza- 
tion factor Z{k) (Figs l(a)-(d)). One sees immediately 
that whatever the adiabaticity, the minimum of E{k) is 
at k = for small A, but at finite k for large A. At first 
glance, nevertheless, nothing unusual seems to happen to 
the ground state energy EqsW at the critical value A c , 
where kcs first becomes non-zero (Fig. 1(e)). In fact, 
the curves in Fig. 1(e) look quite similar to those for 
Holstein polarons. 

However, there is actually a singularity at A c . Plots 
of the dimensionless derivative dEcs(\)/da (Fig. 2(a)), 
the overlap ZqsW between the ground state at finite 
A and the uncoupled ground state (Fig. 2(b)), the mo- 
mentum kcsW for which E{k) is minimized (Fig. 2(c)), 
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FIG. 2: (Color online) (a) Derivative of the GS energy with 
respect to a, (b) Z-factor of the GS, (c) wave vector of the GS, 
and (d) the ratio m D /m* of the bare and effective polaronic 
masses at kas for u p h = 0.5 (here m = l/2t ). Red triangles, 
green rhombi, black squares and blue circles correspond to 
LPBED, MA, DMC, and BDMC methods, respectively. The 
vertical dashed arrow indicates the critical coupling A c . 

and the renormalized effective polaron mass m* (A) = 
{d 2 E{X)/dk 2 )' 1 \ k=kGS (Fig. 2(d)), all show a sharp 
transition at A = X c (uj p h) (see Fig. 2 for = 0.5, 
Fig. 3 for ujph = 3). At this singularity, the pola- 
ronic mass m*(A) diverges, with corresponding jumps in 
the first derivatives dkcsW/dX and dZcsW/dX, and in 
d 2 Eos(\)/dA 2 . However the average number of phonons 
N p h(X) in the polaronic polarization cloud does not di- 
verge at A c (although it is presumably still singular); 
Nph(Xc) < 15 for all values of the adiabaticity param- 
eter ujph checked so far. Note also how A c varies with 
Wj,;, (Fig. 4), initially increasing for small uj p h, but then 
falling to the asymptotic value A c — > 1/2 in the instanta- 
neous phonon limit u)^ oo. This limit can be derived 
analytically (see below). 

We emphasize here the remarkable agreement obtained 
between all 4 methods. The three numerical techniques 
are in principle exact, but all have their practical limi- 
tations, such as the sign problem noted below for QMC 
methods. 

(ii) Discussion: The key new feature of couplings 
M(k,q) like that in Eqs. ((2|) and ((3|), compared to k- 
independent couplings, is that they are non-diagonal in 
site index. Thus phonons cause the bandwidth to fluc- 
tuate, and can by themselves generate hopping between 
sites. The lowest-order process contributing to Egs is 
2nd order in M(k, q); higher corrections come from even 
powers of M. The same applies to the polaron self- 
energy, the polaron mass, quasiparticle renormalization, 
etc. Consider now a pair of vertices, connected by a 
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FIG. 3: (Color online) The same as Fig. [2] but for u v h = 3. 
MA results in (b), (c) are shown as green solid lines. 
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doubly suppressed, because each requires an intersite po- 
laron cloud overlap. Here, however, an electron can hop 
from i — 1— > i — + using only V, to first create 
and then remove a phonon at site i. The associated 
energy is t\ oc M iti+ iAt it i-i/u p h where AU,i-i ~ a 
is the phonon-induced change in the hopping. Since 
the phonon-induced displacement Xi increases one bond 
length while decreasing the other, At^i+iAt^j—i < 0, 
ie., t\ oc —a 2 /uj p h is negative, favoring a minimum in 
E(k) ~ —2^2 cos(2fc) at fc = 7r/2, consistent with our re- 
sults for large A. This simple analysis indicates how the 
transition can occur. Of course, higher order terms must 
also be considered; and a transition like this, signalled 
by the change in kes, is certainly not guaranteed for all 
fc-g-dependent couplings (thus the Edwards model in the 
large A limit also has a dominant iSterm of similar origin; 
but t* 2 > 0, and k GS = for all A 

(c) Finally, consider the limit uj p h — > oo for a fixed 
\,t Q . The phonon propagator tends to its static limit: 
D(q,u) = -2oj ph /(uj 2 h - w 2 ) -> D = -2/u) ph . The 
polaron propagator is then dominated by the 2nd-order 
correction in V, scaling like a 2 /uj p h ~ Ai D (higher order 

/ \ n ~ i 



corrections 



to lowest order in w^ 1 we get from Eq. (j4]) that 
E(k) = -2< cosfc + ^ I75J D^|M(fc,g)| 2 

= — 2i Q cosfc =^ y^[sin(fc + o) — sinfcl 2 (5) 



FIG. 4: (Color online) Phase boundary dividing GSs with 
zero and nonzero momentum. Green squares and red triangles 
refer to MA and LPBED methods. The horizontal dotted line 
refers to the instantaneous limit X c ph =1/2. 



phonon of momentum q; we have 
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Three key new features appear in ((4]): 

(a) it can be of either sign when k ^ k! . This leads to 
a 'sign problem' in any Monte Carlo calculation (indeed, 
for any interaction M(k,q) with non-definite sign); we 
discuss this in the supporting material. The SSH model 
is thus representative of a large class of models in which 
non-diagonal couplings give a sign problem. 

(b) multi-site hopping terms involving phonons gen- 
erate terms in the polaron dispersion of form E(k) = 
Eq — 2t\ cos k — 2^ cos(2fc) — .... Now for g-only depen- 
dent couplings, the nearest-neighbor hopping f* <C t a 



is exponentially suppressed, and t\ 



4At 



18] is 



We see that the dispersion curvature [d 2 E(k)/dk 2 ~\ k _ = 
4t a (l/2 — A) at k — 0. Thus, in the large uj p h limit, the 
effective mass diverges for \ c (u) p h — > oo) = 1/2. Fig. 4 
shows it converges very slowly to this limit. 

This discussion shows, at least for sufficiently large 
u) p h, that there must be a critical coupling strength A c 
at which kes leaves zero. For small to p h the existence 
of a critical point is less clear, because the higher order 
diagrams can have arbitrary sign; but it is what we find 
here for all u: p h studied [l9(. Note, however, that for 
w p d < 0.3, the average number of phonons N p h increases 
significantly, making numerical simulations very difficult. 
The MA method is also questionable in this limit. 

One is tempted to call this T = transition a 'quan- 
tum phase transition'. However this is not correct, be- 
cause any phase transition must involve the cooperative 
behaviour of an infinite set of degrees of freedom; but here 
the number of phonons N p h in the polaronic cloud always 
remains small. Of course with a macroscopic number of 
polarons in the system, we would see non-analyticity in 
bulk properties like dEcsW/dX; but a small number of 
polarons will be invisible in any thermodynamic property. 
Thus we simply assert the existence of a non-analyticity, 
as a function of A, in the polaronic properties. 
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We see that polarons having a coupling to a bosonic 
field depending on both k and q behave in a fundamen- 
tally different way from the standard case with only q- 
dependent coupling. This suggests a large zoology of 
so far unexplored behavior in many physically relevant 
systems. Note how surprisingly different the polaronic 
properties are here. For example, for large A, m*(A) de- 
creases, and Z(kos) remains quite large. We see that 
"standard polaronic behavior" is really just a feature of 
models like the Holstein and Frohlich model. 

Experimental signatures of the new behavior - notably, 
the critical point - will clearly be invisible in any thermo- 
dynamic measurements. However the divergence of the 
effective mass should be easily detectable in transport 
measurements; the polaron mobility fi ~ 1/m* goes to 
zero at the critical point. Thus in any system where the 
charge mobility is carried by the polarons, this critical 
point should be very obvious. It would also be interest- 
ing to do ARPES experiments [20| . where polarons can 
be ejected directly from the insulating state, allowing di- 
rect measurement of E(k) and Z(k). Apart from poly- 
acetylene, various organic semiconductors are known to 
have important non-diagonal coupling to phonons [2lj |. 
as do several dimerized Mott magnetic semiconducting 
oxides [22j; in some of these, the coupling can be varied 
somewhat by pressure. However, any quantitative the- 
ory for such experiments must also include the coupling 
to longitudinal phonons, electron-electron interactions, 
and inter-chain coupling. 
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SUPPORTING MATERIAL 

(i) Quantum Monte Carlo Methods: As discussed 
in the main text, a key feature of the problem addressed 
here is that any Quantum Monte Carlo (QMC) method is 
faced by a sign problem (whereas in the conventional po- 
laron problem, with g-dependent coupling only, the imag- 
inary time polaron Green's function and self-energy are 
sign-definite). Since the sign change can only happen 
when the internal particle line changes momentum, we 
see that non-crossing diagrams are positive-definite and 
only a fraction of the crossing diagrams will have an over- 
all negative sign. This is of course not as severe as the 
sign-alternating series found in the many-body fermionic 
case. 

The Bold Diagrammatic Monte Carlo (BDMC) tech- 
nique [1] is a sign-problem tolerant method for many- 
body problems. It is thus a natural fit to the SSH po- 
laron, and is one of the four techniques used in this letter. 
It improves on the simpler Green's Function Diagram- 
matic Monte CarloQ (G-DMC)and Self-Energy Monte 
Carlo (S-DMC) methods. The former is a Monte Carlo 
sampling of the conventional Feynman diagram expan- 
sion for the Green's function G(k, r) where k and r are 
the momentum and imaginary time respectively. The 
S-DMC method also expands the self-energy E(fc, r) di- 
agrammatically, but contains fewer terms, thus provid- 
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ing faster convergence. The BDMC method further re- 
duces the number of diagrams to be sampled by self- 
consistcntly renormalizing the particle propagator, and 
using this to draw subsequent diagrams. By repeated 
iteration of this procedure, the number of diagrams ac- 
counted for by the BDMC method, grows exponentially 
with simulation time, instead of linearly as in other QMC 
methods. 

In both E-DMC and BDMC the renormalizcd propa- 



gator G'(k, t) is obtained by solving Dyson's equation in 
imaginary frequency £ using a FFT algorithm to convert 
between r- and ^-dependent quantities. Normalization of 
S(fc, r) is enforced using ^ Q T,(k, r) = \ut[2- cos (2k)]. 
Quantities such as the polaron dispersion, quasiparticle 
weight and effective mass can all be obtained either from 
r) or G(k, r). Fig. [5] compares the diagrammatic 
formulation of the various QMC methods. 
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FIG. 5: Comparison between the (a) G-DMC, (b) E-DMC and (c) BDMC methods. Dashed-lined boxes represent MC sampling 
of the appropriate diagram expansion. Some extra restrictions apply in the case of BDMC to avoid double counting. Step 2 
and 3 of BDMC are repeated until desired convergence. Step 3 consists in solving Dyson's equation in imaginary frequency. 



(ii) The Momentum Average Method: The Mo- 
mentum Average (MA) approximation solves the equa- 
tion of motion for the Green's function G(k, uS) — 
(0|cfc(w + ir) — "H) _1 c^|0), where |0) is the vacuum. As 
usual, this solution appears in the form of an infinite hi- 
erarchy of coupled equations of motion. The main idea 
behind MA is to simplify the coefficients multiplying the 
generalized Green's function in these equations so that 
this infinite hierarchy can be solved (quasi) analytically. 
Because it is formulated in terms of equations of motion 
and because there is no truncation at any order, MA 
is a non-perturbative approximation. In diagrammatic 
terms, MA is equivalent [3] to a sum of all diagrams for 
the self-energy expansion, after discarding exponentially 
small contributions (compared to what is kept) from each 
diagram. As a result, MA satisfies multiple spectral 
weight sum rules exactly. 

The guidance on how to approximate the equations 
of motions (which also determines what parts of the dia- 
grams are discarded) comes from the variational meaning 
of MA. Essentially, one makes some assumption about 
the nature of the bosonic cloud, and keeps only terms 
consistent with it For the Holstein model, it is al- 

ready a very good approximation to assume that all the 
bosons in the cloud are at the same site - however, there 
can be arbitrarily many such bosons, and the particle 
can be arbitrarily far from the cloud. Further system- 



atic increases in the variational space lead to systematic 
improvements of the MA approximation [3j. For more 
complicated models, one needs to allow the cloud to be 
more extended. In this work, we use the straightforward 
generalization to our "SSH model" of the 3-site version 
of MA detailed for the Edwards model in Ref. Q . 

(iii) The Limited Phonon Basis Exact Diago- 
nalization method: The Limited Phonon Basis Ex- 
act Diagonalization (LPBED) methodQ is based on the 
standard Lanczos algorithm. It uses the translational 
symmetry associated with periodic boundary conditions, 
requiring that the states have a definite momentum (the 
Hamiltonian is block diagonalized) . Each basis vector 
is a linear superposition with appropriate phases of the 
translational copies (charge carrier and lattice configura- 
tions are together rigidly translated) of a state having the 
electron fixed at a site and phonon quanta located around 
it. The real bottleneck comes from the Hilbert space re- 
quired by the phonon basis, which is unlimited even in a 
finite size lattice. To circumvent this difficulty, LPBED 
takes into account only a finite number of phonon states. 
The states are chosen starting from the observation Q 
that the MA approximation within the Holstein model 
can be recovered by using a restricted basis where all the 
bosons in the cloud are at the same site and arbitrarily far 
from the electron. In the SSH model, where the coupling 
charge-lattice is not local and is related to the hopping, 
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the phonon cloud has to be more extended to describe in 
an appropriate way the relevant physical processes (the 
extension increases as uj ph decreases). In this work we 
use a cluster of 5 neighbouring sites (the cluster can be 
arbitrarily far from the electron). Moreover we include 
an additional pair of phonons located on any pair of lat- 
tice sites. In this way the scattering processes between 
the charge carrier and up to two phonons in q-space are 
exactly treated, so that the weak coupling limit is op- 
timally described. The main advantage of the LPBED 
method is that it is possible to calculate not only the self 
energy of the quasiparticle but any correlation function. 
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